Horizon-absorption effects in coalescing black-hole binaries: An effective-one-body 

study of the non-spinning case. 
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We study the horizon absorption of gravitational waves in coalescing, circularized, nonspinning 
black hole binaries. The horizon absorbed fluxes of a binary with a large mass ratio (q — fOOO) 
obtained by numerical perturbative simulations are compared with an analytical, effective-one-body 
(EOB) resummed expression recently proposed. The perturbative method employs an analytical, 
linear in the mass ratio, effective-one-body (EOB) resummed radiation reaction, and the Regge- 
Wheeler-Zerilli (RWZ) formalism for wave extraction. Hyperboloidal (transmitting) layers are em- 
ployed for the numerical solution of the RWZ equations to accurately compute horizon fluxes up to 
the late plunge phase. The horizon fluxes from perturbative simulations and the EOB-resummcd 
expression agree at the level of a few percent down to the late plunge. An upgrade of the EOB 
model for nonspinning binaries that includes horizon absorption of angular momentum as an ad- 
ditional term in the resummed radiation reaction is then discussed. The effect of this term on the 
waveform phasing for binaries with mass ratios spanning 1 to 1000 is investigated. We confirm that 
for comparable and intermediate-mass-ratio binaries horizon absorbtion is practically negligible for 
detection with advanced LIGO and the Einstein Telescope (faithfulness > 0.997). 
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I. INTRODUCTION 

The dynamics of the quasi-circular inspiral of coalesc- 
ing binary black hole (BBH) systems is driven by the loss 
of mechanical angular momentum through gravitational 
radiation. The total loss of angular momentum consists 
of two contributions: the one due to radiation emitted to 
future null infinity (T^), and the one due to radiation 
absorbed by the black- hole horizons (J 7 ^)- Typically the 
former dominates over the latter, i.e. T$ > . For 
example, the leading order contribution to horizon ab- 
sorption for a nonspinning binary is a 4PN contribution 
of the form [1-3] 



y 



H 



x 4 (1 - 4i/ + 2v 2 ) [1 + <${y)x + 0(x 2 )] . (1) 
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Above x = (AIQ) 2 / 3 is the post-Newtonian orbital pa- 
rameter, f2 is the orbital frequency, M = Ma + Mb is 
the total mass of the system, with Ma,b the masses of the 
individual black-holes, v — MaMb/M 2 is the symmetric 
mass ratio, and = v 2 32/5 x 7 / 2 is the Newtonian con- 
tribution to the asymptotic flux. The explicit expression 
of Ci (v) follows from the state-of-the-art lPN-accurate 
result of Taylor and Poisson [2]. In the presence of spin, 
a more complicated formula holds [2] , with the contribu- 
tion of absorption entering already as a 2.5PN effect. In 
practice, horizon absorption is a negligible effect when: 
(i) the separation between the two objects is large: (ii) 
the two objects have comparable masses (y ~ 1/4); (hi) 
the spins are small. 

Leading-order calculations by Alvi [1] (improved to 
1PN fractional accuracy by Taylor and Poisson [2]) es- 
timate the effect of horizon flows on the number of grav- 



itational wave (GW) cycles to be no more than 10% of 
a cycle for comparable-mass (q — Mb /Ma = 4) bina- 
ries with maximally spinning black holes by the time of 
merger (see Table IV of Rcf. [1]). In the nonspinning 
case absorption effects seem negligible with accumulated 
dephasings that are smaller than 1% of a cycle. 

The analysis of [1] is, however, inaccurate during the 
late inspiral and plunge (1/6 < x < 1/3). In this regime, 
absorption effects may be relevant for GW detection due 
to relativistic corrections, if the mass ratio or the individ- 
ual spins are sufficiently high. The potential importance 
of absorption effects during the late plunge of spinning 
binaries was also pointed out by Price and Whelan [4] 
using the close-limit approximation. 

To meaningfully ascertain the importance of energy 
and angular momentum flows in or out of the black holes 
(depending on the orientation of the spin with respect to 
angular momemtum) during the late inspiral and plunge, 
one needs numerical relativity (NR) simulations. The 
growth rate of the irreducible mass and angular momen- 
tum of the black hole horizons in a NR simulation of 
nearly-extremal spinning black hole binary [5] has been 
compared to Alvi's analytical prediction. A remarkable 
numerical agreement between the two was found up to 
x < V0.5 ~ 0.71, while significant deviations from nu- 
merical data were observed for larger values of x. This re- 
sult suggests that horizon-absorption effects should be in- 
corporated in the analytical modeling of coalescing black 
hole binaries. 

To bridge the gap between the leading-order estimate 
of Alvi valid during the early inspiral [1] and the qualita- 
tive understanding of Price and Whelan valid during the 
late plunge [4] one needs an analytic description of the 
absorbed fluxes that incorporates high-order PN correc- 
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tions and that is not limited to the slow-velocity, weak- 
field regime. Focusing on nonspinning binaries, Ref. [3] 
adapted the resummation procedure of the asymptotic 
energy flux of Ref. [6] to the energy flux absorbed by the 
two black holes, so to consistently incorporate it within 
the effective-one-body (EOB) [7-9] description of the dy- 
namics of black hole binaries. The final outcome of that 
study is an analytical expression of the absorbed energy 
flux, written in a specific factorized and resummed form, 
that is well-behaved (contrary to a standard, PN ex- 
pansion) also in the strong-field-fast- velocity regime (no- 
tably, also along the EOB-dcfincd sequence of unstable 
circular orbits). The input for the resummation proce- 
dure is given by state-of-the art PN-expanded results 
for the horizon flux: the 1PN accurate expressions of 
Taylor and Poisson [2] (valid for any mass ratio), and 
the leading-order results of Poisson and Sasaki [10] in 
the test-mass (y = 0) limit. In addition, this analyti- 
cal knowledge was further improved by adding higher- 
order (effective) PN coefficients extracted from the ab- 
sorbed fluxes from circular orbits computed numerically 
in the test-mass limit. Finally, v = and v ^ (scmi)- 
analytical results were hybridized to get improved accu- 
racy for any mass ratio. 

In this paper we study the effect of horizon-absorption 
on the phasing of circularized, coalescing black-hole bina- 
ries up to merger. We do this by using the EOB descrip- 
tion of the binary dynamics and radiation [7-9] . The ra- 
diation reaction is improved by an additional term, , 
that takes into account the loss of mechanical angular 
momentum due to horizon absorption. As a first cut of 
the problem, we consider here nonspinning binaries only, 
where the effects are weaker than when the BHs are spin- 
ning . 

First of all, we focus on the "large-mass-ratio" limit 
(e.g., v = 10~ 3 ) and we check the consistency of the 
{v = 0) analytical expression of (x; 0) proposed in [3] 
against the absorbed GW flux obtained numerically using 
a Rcgge-Wheclcr-Zcrilli (RWZ) perturbative treatment. 
This gives further confirmation of the reliability of the re- 
summation and hybridization procedure of the absorbed 
flux introduced in [3]. Then we perform a comprehen- 
sive EOB study to investigate the effect of F^(x; v) on 
the phasing up to merger with 10~ 3 < v < 1/4. Note 
that NR simulations for mass ratios q = 100 are cur- 
rently doable [18-20], though they are challenging, do 
not yet provide sufficiently long waveforms, and it does 
not seem practical to cover the parameter space densely 
with full numerical relativity simulations only. There- 
fore, the EOB model is of fundamental importance to 



1 Note that the EOB approach can account consistently for (arbi- 
trary) spins [9, 11-14]. Black hole absorption has already been 
included in EOB-based evolutions of extreme-mass-ratio (EMR) 
inspirals around a Kerr black hole, though only in its standard 
Taylor-expanded form [15, 16]. An improved treatment of this 
problem is currently under development [17]. 



investigate the so-called intermcdiate-mass-ratio (IMR) 
regime [21-25]. 

The RWZ time-domain perturbative method employed 
in this work to obtain large-mass-ratio waveforms is de- 
scribed in detail in [26-30]. We solve the RWZ equa- 
tions for a binary system made of a point-particle on a 
Schwarzschild background and subject to leading-order, 
0(v), EOB-rcsummed analytical radiation reaction. The 
main technical improvement introduced here is the de- 
velopment of smooth, transmitting (hypcrboloidal) lay- 
ers [31] attached to a compact domain of Schwarzschild 
spacctime in standard coordinates to include both future 
null infinity, , and the black-hole horizon, H, in the 
computation. With this method, the absorbed and radi- 
ated fluxes can be computed very accurately. Also, the 
finite differencing order has been improved to 8th order 
accurate operators. These technical developments lead to 
such an efficient code that tail decay rates for the late- 
time of the gravitational waveform emitted by inspiraling 
point particles can be computed accurately (this was not 
possible previously using standard methods). 

This paper is organized as follows. In Sec. II we review 
the results of Ref. [3] that are relevant for this work and 
we give the explicit expression for J 7 ^ . In Sec. Ill we 
discuss the construction of transmitting layers and their 
advantages in improving the accuracy of the numerical 
solution of the RWZ equation. In Appendix A we also 
demonstrate that the layer technique helps solving a pre- 
viously difficult problem of obtaining accurate power law 
tails for inspiraling particles. In Sec. IV we present the 
RWZ calculation of the absorbed waveforms and flux and 
the consistency check of J-^ (x ,0). The main results of 
the paper are collected in Sec. IV C, where we investigate 
the influence of J-^ on the phasing up to merger. Con- 
cluding remarks are gathered in Sec. V. We use units 
with G = c = 1. 



II. EOB DYNAMICS AND WAVEFORM: 
INCLUDING HORIZON ABSORPTION 

In this Section we review the main elements of the EOB 
approach and we recall the results of [3] that are needed 
to compute Th- The EOB analytical description of the 
dynamics of a circularized binary essentially relies on two 
building blocks: the resummed EOB Hamiltonian H^ob, 
which describes conservative effects, and the resummed 
mechanical angular momentum loss which describes 
nonconservative effects due to loss of GW energy (radia- 
tion reaction) 2 . The EOB Hamiltonian depends only on 
the relative position and momenta of the binary system. 



2 An additional radiation reaction term, T r , is present due to linear 
momentum loss through GWs, but, for circularized binaries, is 
typically not included because it remains small up to the late 
plunge. 
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For nonspinning binaries it has the structure 



H 



EOB I 



{r,p r „p v ) = M^l + 2v{H eS -l), (2) 



where 



H, 



eS 



i 



A(r) 



Here Z3 = 2^(4 — 3v) and we use rescaled dimension- 
less variables, namely r = tabc 2 /(GM), where tab = 
\ta — the relative separation between the two bod- 
ies, and p v = P v /(GMaMb), the angular momentum. 
In Eq. (3), p rt is the radial momentum canonically con- 
jugate to a EOB-defined tortoise coordinate, r*, that re- 
duces to the usual tortoise coordinate when v = 0. The 
function A(r) is the basic radial potential that, follow- 
ing Rcf. [32], depends on two EOB flexibility parame- 
ters (05, ae) that take into account effective 4PN and 
5PN corrections to the conservative dynamics. For co- 
alescing black-hole binaries, an excellent phasing agree- 
ment between NR and EOB waveforms can be reached 
over banana-like regions in the (05,06) plane. Follow- 
ing Ref. [32], we fix the EOB parameters as 05 = —6.37 
and ae = 50 which lie within the extended region that 
yields a good fit with NR data for q = 1, 2, and 4. A re- 
cent study [33] comparing an (05, 06)-paramctrizcd EOB 
model with NR simulations for q = 1, 2, 3, 4, and 6 (and 
more accurate than those used in Ref. [32]) pointed out 
that the "best fitting" region in the (05,05) plane actu- 
ally depends on v (see Fig. 5 in [33]). Since our goal here 
is to highlight only the effect of on the dynamics, we 
neglect this further ^-dependence on (05,05). The anal- 
ysis of the ^-dependence of (05, a^) in the calibration of 
the EOB model of Ref. [32] , in the presence of black- hole 
absorption and with better numerical data, is postponed 
to future work. 

The radiation reaction force, J- v , drives the angular 
momentum loss during evolution. The Hamilton equa- 
tion for p v reads 



dp v 
~dT 



(4) 



where T v = T^jv. The mechanical angular momentum 
loss is typically written as 
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i/r*n 5 /(«„; v) 



(5) 



Here, Q — dip/dt is the orbital frequency, with ip the 



orbital phase, v v 



, 51 is the azimuthal velocity, and 



r u = rip 1 / 3 , where ip is a ^-dependent correction factor 
that is necessary to formally preserve Kepler's law during 
the plunge [34]. The function f(x; v) is the reduced flux 
function that is defined, for a circularized binary, as the 
ratio between the total energy flux and the I = m = 2 



asymptotic energy flux. In our case the reduced flux func- 
tion is given by the sum of an asymptotic and a horizon 
contribution as 

f{x;u)=f' f {x;v)+f B {x;u), (6) 

where each term is given by 

f^ H \x;v) = F^f/F^. (7 ) 

Here, F^'^ are the total asymptotic {^) and horizon 
(H) energy fluxes for circular orbits summed up to mul- 
tipole I = £ max , while F^ 2 = (32/5)^ 2 x 5 is the Newto- 
nian quadrupolar (asymptotic) energy flux. In the EOB 
model one uses suitably factorized expressions for the 
multipolar fluxes F) ' to resum and improve them 
with respect to standard PN-expanded expressions in the 
strong-field, fast-velocity regime (1/6 < x < 1/3). The 
resummation of the asymptotic waveform and fluxes was 
discussed in Ref. [6] and has been used in many works 
since then. We use it here at the 3 +2 PN accuracy 3 and 
we fix i max = 8. 

The horizon flux is written as the sum (up to £ max = 8) 



(8) 



1=1 m=l 



where the partial multipolar fluxes have the following 
factorized structure [3] 



Seff(z;^) [Plm{x] is))' 



(9) 

Here, e = ir(£ + m) = 0, 1 is the parity of the considered 
multipole, is a source factor, with = H e s or 
iSgg = \fxp v according to the parity of the multipole, 
and the pf m (x; v) are the residual amplitude corrections 
to the horizon waveform. Only /of^>(x; v) is known ana- 
lytically at 1PN accuracy [3]. It reads 



To improve our knowledge of the strong-field behavior 
of the P22 functions, Ref. [3] computed numerically the 
Pim 1 ™ 1 functions for a test particle moving on (stable and 
unstable) circular orbits on a Schwarzschild background. 
For each multipole, it was possible to fit the numerically 
computed p^ um accurately via a suitable rational func- 
tion of the form 



Hm _ 1 + n{ m x + n e 2 m x 2 + nf"x 3 + nfV 



1 + d\ m x + d^x 2 



(11) 



3 The 3PN-accurate ^-dependent terms are augmented by the 4PN 
and 5PN accurate v = corrections for all multipolcs. 
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TABLE I: Coefficients of our hybrid l +3 PN-accurate 
Plmip'i v ) functions as given by Eq. (13). 



1 


m 


lm 
c l 


lm 
C2 


lm 
c 3 


lm 

c 4 


2 


2 




4.78752 


26.760136 


43.861478 




2 


1 


0.58121 


1.01059 


7.955729 


1.650228 


3 


3 


1.13649 


3.84104 


45.696716 


27.55066 


3 


2 


0.83711 


1.39699 


23.638062 


-1.491898 


3 


1 


1.61064 


2.97176 


10.045280 


15.146875 


4 


4 


1.15290 


4.59627 


55.268737 


13.255971 


4 


3 


0.96063 


1.45472 


43.480636 


-35.225828 


4 


2 


1.43458 


2.43232 


21.927986 


10.419841 


4 


1 


0.90588 


1.17477 


5.126480 


4.022307 



where nf 171 and d[ m are free fitting parameters 4 . By 
Taylor-expanding Eq. (11) in powers of x one obtains 
the following representation of the pf m functions in the 
v = limit 

P? m (x;0)=T N [p?«*(x)}, (12) 

where N indicates the maximum power of the expan- 
sion. For the I = m = 2 mode, Ref. [3] pointed out 
that setting N = 4 (i.e., 4PN accuracy) is sufficient to 
yield an accurate representation of the p^ um up to and 
below the last-stable-orbit (LSO) at r — 6, with rela- 
tively small differences around the light-ring (see Fig. 3 
of [3]). We have verified that this remains true also for 
the other multipoles, so that we shall assume 4PN ac- 
curacy in Eq. (12) from now on. Following Ref. [3], we 
hybridize the ^-dependent 1PN information of Eq. (10) 
with the 4PN expansion of Eq. (12). Such hybridization 
procedure, that is conceptually analogous to what has 
been done in Ref. [6] for the corresponding asymptotic 
residual amplitude corrections, is justified in view of the 
following two results of Ref. [3] : (i) the dependence on v 
of the 1PN coefficient in Eq. (10) is mild; (ii) the fit of the 
numerical data proved to be robust enough so that the 
coefficients of the PN expansion can be taken as reliable 
estimates for the actual (yet un-calculated) PN coeffi- 
cients. In practice we use the following 4PN expression 
for the pf m (x; v) 

H f \ 1 i tm m i lm 2 i Ira 3 , 4 /i o\ 

Pe m {x; v) = 1 + cj x + c 2 x + c 3 x + c 4 x . (13) 

The values of the coefficients cf Tn , i = 1, . . . , 4 are listed 
in Table I, where in fact only c 22 is given analytically as 
a function of is, while the other coefficients are computed 
from the test-mass nf 71 and df m coefficients extracted 



4 Note that for the £ = m = 2 mode the fit was done imposing 
the constraint that the 1PN coefficient is equal to 1, because 
pg(x; 0) = l + x + 0{x 2 ). 



from the fit. We shall use them in the following as an 
effective representation of the actual test-mass informa- 
tion, although the hope is that it will be soon possible to 
replace them with terms from a PN calculation. 

In Table I we list all PN coefficients up to I = 4. It 
seems enough to include only the quadrupolar contribu- 
tions P21 and P22 in f (x; v), since, as we show in Sec. IV 
below, the effect of multipoles with I > 3 on the horizon- 
absorbed angular-momentum flux is practically negligi- 
ble already in small-mass-ratio coalescence events. [Note 
that the independence of the leading-order prefactor to 
the multipolar horizon flux, is fully known only for 

the quadrupole modes [2]]. 

Using Eqs. (6), (9) and (13) one defines an EOB dy- 
namics that takes into account horizon absorption. From 
this dynamics one then computes the (asymptotic) EOB 
multipolar waveform that has the well known factorized 
structure 

hi m = h[ N Js^h\fM m ) l h^ C , (14) 

where h^ e ^ is the Newtonian waveform, /i tal1 = Tf m e I(5fm 
is the tail factor as defined in Ref. [6], pi, n is the re- 
summed modulus correction and is a next-to-quasi- 
circular correction. For each multipole (£, m) these NQC 
corrections depend on 4 parameters, a\ m , i = 1,...,4 
(two for amplitude corrections and two for a phase 
correction) that have to be determined with an itera- 
tive procedure to match the EOB waveform to the NR 
waveform around merger. The NQC correction to the 
amplitude depending on (af 11 ,^ 171 ) is the same as in 
Refs. [29, 32]; the NQC correction to the phase depending 
on (o| ro , af 1 ) is implemented as per Eq. (22) of Ref. [33], 
that proved more robust than the analogous expression 
used in Eq. (12) of [29] to complete the EOB waveform 
in the extreme mass-ratio limit. The a\ m parameters are 
determined as in [29] by imposing that the slope of the 
EOB waveform amplitude and frequency agree with the 
NR ones at the peak of the EOB orbital frequency fi. 
Note that, consistently with the findings of [29] and dif- 
ferently from previous work [32, 33], we do not impose 
that the peak of |/i22 1 occurs at the same time as the peak 
of Cl. On the contrary, we allow I/122I to have a nonzero 
slope there that coincides with the slope of the NR wave- 
form modulus I/I22 I at a NR time that occurs slightly 
after the time corresponding to max |/i22|- This NR-data 
extraction point is suitably chosen consistently with the 
test-mass results [35]. To obtain the coefficients af m for 
any value of v, we fit with cubic polynomials in v the 
NR points extracted from both the waveforms computed 
for us by D. Pollncy and C. Rcisswig using the Llama 
code [35-37], for mass ratios q = 1,2,3,4, and the per- 
turbative data of [29, 30]. As a last step we match to the 
EOB inspiral-plus-plunge waveform, Eq. (14), a super- 
position of Kerr black-hole quasi-normal-modes (QNMs) 
over a matching "comb" [27]. We use in general five 
QNMs; note, however, that for v = three QNMs are 
sufficient to obtain good agreement between EOB and 
RWZ waveforms [29]. 
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III. TRANSMITTING LAYERS FOR THE 
REGGE-WHEELER-ZERILLI EQUATIONS 

In this Section wc describe the transmitting layers 
adopted here to solve the RWZ equations and to extract 
the GW fluxes at the horizon and at null infinity. The 
method builds on previous work [30, 31, 38] and extends 
the hypcrboloidal layer technique. We also present, as a 
test of the implementation, the horizon absorbed fluxes 
from geodesic circular motion, and, in Appendix A, we 
report tail computation at with our new infrastructure 
based on transmitting layers. 



in this layer is then not hypcrboloidal but horizon pene- 
trating. Therefore, we occasionally use the more general 
term transmitting for this new layer. The name is chosen 
to distinguish it from absorbing layers commonly used in 
applied mathematics [40]. Instead of artificially absorb- 
ing the outgoing signal, our layers transmit it to infinity, 
hence the name. 

The method consists of a spatial coordinate compact- 
ification and a time transformation as described below. 



1. Spatial compactification. 



A. Smooth transmitting layers 

We use the Schwarzschild time coordinate t and the 
tortoise coordinate r* in the bulk for describing the inspi- 
ralling particle using the standard EOB formalism. The 
tortoise coordinate 



r, = r + 2Mlog(r-2M), 



(15) 



is constructed such that the event horizon r = 2M is 
at infinite coordinate distance. From a numerical point 
of view, the main effect of the tortoise coordinate is to 
push away the coordinate singularity at the bifurcation 
sphere in Schwarzschild coordinates. The computational 
domain is then truncated at some negative value for r» 
and ingoing boundary conditions are applied. 

There are two problems with this common approach. 
First, the artificial truncation of the computational do- 
main leads to artificial boundary conditions. This prob- 
lem is not as important in the negative r* direction as 
in the positive one, because the potential falls off expo- 
nentially in the tortoise coordinate towards the horizon 
whereas only polynomially towards spatial infinity. Nev- 
ertheless, the imposition of such artificial boundary con- 
ditions can still complicate the implementation of higher 
order discretization methods. Second, the computation 
of absorbed fluxes by the black hole is performed at finite 
radius. To avoid contamination of the horizon flux com- 
putation by the artificial boundary conditions, a large 
grid in the negative direction needs to be chosen (see, 
for example, [39]). This practice leads to a waste of com- 
putational resources. 

A resolution to these problems is to change the 
coordinates near the horizon and in the asymptotic 
domain ("near infinity"), while keeping the standard 
Schwarzschild coordinates in the bulk. In our previous 
studies [29, 30] we applied hypcrboloidal scri-fixing in a 
layer [31, 38] to solve these problems near infinity. In its 
original form, such a hypcrboloidal layer is attached in 
the positive radial direction so that the outer boundary 
corresponds to future null infinity. Since we are using the 
tortoise coordinate r% , a similar layer can be attached also 
in the negative r + direction so that the inner boundary 
corresponds to the black hole horizon. The time foliation 



Consider a finite domain T> in the tortoise coordinate 
r* given by V = [— i?_,i?+] where R± £ M + . In this 
finite domain, we use coordinates (t, r*). Wc introduce a 
compactifying coordinate 5 p to calculate the solution to 
the RWZ equations numerically on the unbounded do- 
mains (— oo,— R-) and (i?+,oo). The compactification 
is such that the infinities are mapped to a finite p, and 
at the interfaces R± the coordinates p and agree. 

A convenient way to write such a compactification is 



(16) 



where Q(p) is a suitable function of p. It is unity in 
the bulk domain, fip = 1, implying p = r* on T>. For 
compactification, £1 must vanish at a finite p location, 
which then corresponds to infinity with respect to (see 
[30, 31] for details). The transformation therefore is de- 
generate at the zero set of f2. Its Jacobian reads 



J 



dp n 2 
~ n - p w '■ 



(17) 



where the prime indicates d/ dp. A simple prescription 
for f2 to compactify both directions could be 



n = i 



\p\ - R± 

S±-R± 



Q(\p\-R±) 



(18) 



For p < we use the plus sign, for p > we use the minus 
sign in the above formula. The transformation (16) with 
(18) maps the unbounded domain — oo < r* < +oo to 
the bounded domain — S- < p < 5+ such that p = r» on 
V = [-R_,R+\ where S± > R±. 

The choice of f2 in (18) leads to a coordinate transfor- 
mation that is C 4 at the interfaces. Our numerical ex- 
periments showed that this degree of smoothness was not 
sufficient for the accurate computation of late-time tail 
decay rates of the waveform reported in Appendix A. Nu- 
merical studies of hypcrboloidal compactification using 



For notational continuity with previous work we use the same 
symbol to address both the compactifying coordinate and the 
residual amplitude corrections pi m to the EOB waveform. 



6 





h6- 






0.8 






0.6 




j 


0.4 






0.2 





-20 



-10 



10 



20 



FIG. 1: The function fi for the two choices (18) and (19). 
The dashed vertical lines (red online) indicate the interfaces 
at R± = 12. Infinity corresponds to S± = 20. The dashed 
(black online) curve denotes the C 4 transition (18), and the 
solid (black online) curve the smooth transition (19). The 
numerical results obtained later in the text are obtained with 
the smooth transition and with this choice of parameters. 



RWZ equations previously showed that a smooth (C°°) 
transition leads to higher accuracy [41]. Such a smooth 
transition function can be given as 



St ■= 2 + \ tanh 



where we have defined 



— ( tan x — 

it V tana; 



7T \p\ 



2S±-R±) 



The free parameter q determines the point pi/% at which 
/t(Pi/2) = 1/2 and s determines the slope of /t at p 1 ^ 2 
[42, 43]. We set 



Q = 1 



H 
s± 



f T Q(\p\-R±). 



(19) 



The two choices for f2 have been plotted in Fig. 1. 
For the main numerical results in this paper we use the 
smooth compactification of Eq. (19). 



2. Time transformation. 



Such time transformations can be written in the follow- 
ing form 



t±h(r*) 



(20) 



where the function h is called the height function and 
depends on the tortoise coordinate only. 

Near the black hole horizon, and near null infinity, 
gravitational waves propagate predominantly in one di- 
rection along null rays. Near the black hole most waves 
are absorbed, near infinity most waves escape. Corre- 
spondingly, near the black hole we require invariance of 
ingoing null rays in local coordinates, whereas near in- 
finity we require invariance of outgoing null rays. The 
sign in Eq. (20) depends therefore on the sign of r*. The 
invariance of the null direction in local compactifying co- 
ordinates translates into 

t ± = t ± p. 



With Eq. (20) we get 

r* = p + h 
or by defining H := dh(r*) / dr* 

H =1 — J, 



(21) 



This relation between the differential time transforma- 
tion H and the differential spatial compactification J 
solves the resolution problem of compactification in hy- 
perbolic equations. 

We emphasize that, even though the inner transmit- 
ting layer changes the time foliation, we do not mod- 
ify the particle trajectory consistently when solving the 
RWZ equation. In principle, the particle motion should 
be expressed in the local coordinates of the inner layer. 
In practice, however, this seems unnecessary when the 
layer is attached at a sufficiently small negative value 
of = — i?_ < 0. We find that after the particle has 
crossed the light ring at 3M, thereby triggering the QNM 
ringdown, its subsequent trajectory does not influence 
the waveform. Choosing i£_ = 12 allows us to leave the 
description of the particle untouched. Once the particle 
enters the layer, we smoothly switch off the RWZ source 
to avoid unphysical features in the ringdown waveform 
(see Fig. 16 of [30]). 



It is well known that spatial compactification alone 
leads to resolution problems for hyperbolic equations [44] . 
The loss of resolution near infinity, however, can be 
avoided for essentially outgoing solutions by combining 
the spatial compactification with a suitable time trans- 
formation [31]. The details of this transformation depend 
on the background spacctime, but the essential idea is to 
keep the outgoing null direction invariant in local com- 
pactifying coordinates [30]. 

A suitable time transformation for numerical computa- 
tions keeps the background metric invariant of the time 
coordinate by respecting the timelike Killing field [38]. 



B. Horizon fluxes for circular orbits 

As a test of the accuracy of our new numerical setup, 
and in particular, of the inner layer, we consider a point- 
particle moving on circular orbits of a Schwarzschild 
black hole and we compute the horizon fluxes. The 
treatment of the distributional <5-function describing the 
point-particle source as a finite-size, narrow Gaussian 
is the same as previous works [26-30]. Given a se- 
lected sample of stable and unstable orbits of radius r 
(3.1 < r < 7.9 spaced by Ar = 0.1), the RWZ waveform 
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FIG. 2: (color online) Testing the accuracy of the updated 
time-domain RWZ code for a particle along a sequence of sta- 
ble and unstable circular orbits. We plot the fractional dif- 
ference in the horizon fluxes computed with the time-domain 
RWZ code using transmitting layers and S. Akcay's frequency- 
domain code [3, 45]. 



(H,e) 



and its time derivat 



lve, 



at the horizon location ^ 

$J m ' , the fluxes of energy and angular momentum ab- 
sorbed by the black hole are given by [39] 



FIG. 3: (color online) Comparing horizon and null infinity 
quadrupolar (I = m = 2) RWZ waveforms for a coalesc- 
ing binary with mass ratio v = 10 -3 . The horizontal axis 
corresponds to horizon anticipated time u = u + = r + H 
for the horizon waveform and to null infinity retarded time, 
u = u~ — t — S for the asymptotic waveform. The leftmost 
(dash-dotted) vertical line marks the (dynamical) time when 
the particle crosses the LSO, while the rightmost (dashed) ver- 
tical line corresponds instead to the light-ring crossing. The 
horizon waveform (red online) becomes unreliable around the 
light-ring crossing {u/M > 4300). See text for discussion. 



IV. HORIZON ABSORPTION IN THE 
LARGE-MASS-RATIO LIMIT 



E, 



J, 



H 



j_ — (i + 2 )! (g<6) 

1Qn U-2)V lm 1 

1=2 m=0 e=0 v ' 



(22) 



1=1 m=l e=0 



(1 + 2)! ^ 



(23) 



In Fig. 2 we show the fractional difference (plotted ver- 
sus x = 1/r = (A/0) 2 / 3 ) between the the energy flux 
E H computed with our code (labeled by "BNZ") and 
the same quantity obtained by S. Akcay using his fre- 
quency domain code [45] , and presented for the first time 
in Ref. [3] (labeled by «NA"), i.e., (£; BNZ - E NA )/E NA . 
The solid (red online) curve in the plot refers to the total 
flux summed up to f max = 8, while the dashed one to 
the £ = m = 2 dominant quadrupole mode only. The 
frequency domain computation of horizon fluxes using 
the code of Ref. [45] have fractional uncertainty of order 
10~ 10 or smaller for strong-field orbits (say r < 10). Fig- 
ure 2 highlights how the fractional difference between the 
fluxes obtained with the two methods is on the order of 
10" 3 . 



A. Perturbative, time-domain computation 

In this section we compute the horizon-absorbed GW 
fluxes in a large-mass-ratio BBH coalescence using the 
perturbative method discussed extensively in previous 
works [26-30]. The computations allow us to test the 
reliability of the EOB-resummed fluxes given by Eq. (9). 

In the large-mass-ratio limit the EOB Hamiltonian 
tends to the Schwarzschild one, and higher-order correc- 
tions in the analytical radiation reaction are neglected. 
The radiation-reaction term is then given by 



T 



H 
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with v v = rfi 



0) + f H (v v ;0) , 
(24) 

Here, y 0) is computed as in Ref. [6] 
in the v = limit but retaining all terms up to 5PN 
fractional accuracy in the pf m 's computed in Ref. [46] 
(see also Ref. [47] for the 14PN accurate calculation). 
We work here with the mass ratio 6 v = 10 -3 . Previous 



6 Note that in the test-mass limit, M^/Mg 1 we can identify 
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studies [28, 30] indicated that, in this case, the method 
gives a fractional agreement between the 5PN-accurate 
mechanical angular momentum loss and the actual an- 
gular momentum flux computed from the RWZ master 
function of order 10 -3 even beyond the LSO (see Fig. 14 
of [30].) The RWZ master function is extracted numeri- 
cally using the method of Sec. III. Neglecting horizon ab- 
sorption in the dynamics {f H {v v ; 0) = 0, in Eq. (24)), we 
reproduce the relative dynamics of previous works [28- 
30]. The initial relative separation is — 7 and the 
relative dynamics is started with the usual post-circular 
initial conditions [8, 26]. 

Figure 3 focuses on the £ = m = 2 mode and illustrates 
the relative importance of the horizon waveform ^> 2 ^ '°' 

compared to the asymptotic waveform 4*22 ■ The fig- 
ure shows on the same panel the real part of the wave- 
forms together with their amplitudes. In the strong-field 
regime under consideration, r < 7, the horizon waveform 
is smaller (~ 16 times during inspiral) than the asymp- 
totic waveform but not negligible (roughly comparable to 
some asymptotic subdominant multipoles). Notably, one 



finds that |* 



22 



is always larger than |^ 44 



(^,o). 



The ra- 



tio between the two varies between 1.5 at the beginning 
of the inspiral up to 2 at LSO crossing. 

The amplitude of the horizon waveform grows during 
the late plunge and reaches about 0.1 just before the 
light-ring crossing, u/M sa 4300. It then increases by a 
factor ~ 7 over a temporal interval ~ 15, developing a 
"spike" that is twice as large as the corresponding value 
of the asymptotic amplitude. After this transient, the 
ringdown asymptotic and horizon waveforms are consis- 
tent. 

The presence of a spike in the horizon waveform is due 
to our representation of the point-particle source as a 
narrow (cr -C 1) Gaussian. The RWZ function is (in the 
o — y limit) discontinuous at r» = R(t) and its spatial 
derivative is singular. Since we have not implemented 
a sophisticated regularization of the source (see in this 
respect Refs. [48-51]), there is a spatial (smoothed) sin- 
gularity on the RWZ computational grid at the particle 
location. After the particle has crossed the light ring, the 
singularity is advected to the horizon. The presence of 
such a discontinuity in the RWZ function and the corre- 
sponding singularity in the energy flux (also observed in 
the analytical treatment of an extreme-mass-ratio plunge 
by Hamerly and Chen [52]), makes our numerical repre- 
sentation of the particle ill-suited for a detailed study 
of horizon absorption during the last moments of the 
merger. We have, however, verified that the effect is 
localized around the location of the particle and its influ- 
ence is reduced for smaller values of a. In this work, we 
use the RWZ horizon waveform (and flux) only before the 
light-ring crossing, say u/M ~ 4300, so that our results 
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the inverse mass ratio l/q = Ma/Mb with the symmetric mass 
ratio v = M A M B /(M A + M B ) 2 . 
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FIG. 4: (color online). Top panel: comparison between RWZ 
horizon and asymptotic angular momentum fluxes for mass 
ratio v = 1CP 3 from Eq. (23) with £ max = 8. Bottom panel: 
the £ = 2 modes contribute to more than the 98% of the total 
absorbed flux up to LSO crossing (vertical dashed line). 



are not affected by the absorption of the particle by the 
horizon. 

We display in Fig. 4 the horizon-absorbed angular 
momentum flux jv 2 computed from Eq. (23) with 
4iax = 8. The top panel contrasts asymptotic fluxes (ei- 
ther summed up to £ max = 8 or just £ = m = 2), with the 
horizon fluxes, highlighting that the latter are typically 
10~ 3 times smaller. The bottom panel of the figure shows 
the ratio between the total quadrupolc horizon flux (i.e., 
J 21 + J 22) and the total horizon flux jE _ g s, which 
indicates that the quadrupolc mode accounts for more 
than the 98% of the absorption up to the LSO crossing 
(dash-dotted vertical line in the plot). 



B. The EOB-resummed horizon flux 

We compare the horizon absorbed angular momentum 
flux computed from the RWZ waveform, Eq. (23), with 
the EOB-dcfincd mechanical angular momentum loss due 
to horizon absorption, Eq. (24). In this section, the dy- 
namics is computed including only T£ ; the effect of 
is explored in the next section. Figure 5 shows the dom- 
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FIG. 5: (color online) Comparison between EOB resummed angular momentum flux and the RWZ one for the quadrupole 
(£ — 2) modes: m = 1 (left panel) and m = 2 (right panel). The dash-dotted vertical line line marks the LSO crossing. 
The EOB-resummed (horizon) mechanical angular momentum loss is perfectly consistent with the horizon flux computed from 
GWs. By contrast, the lPN-accurate expressions, Eqs. (25)-(26), underestimate horizon absorption by more than a factor 2. 



inant quadrupole i = 2 fluxes for m = 1 (left panel) 
and m = 2 (right panel). The mechanical losses —T^/v 
computed with various approximations (non-solid lines) 
are contrasted with J^ m /v 2 (solid lines) The vertical 
dash-dotted line marks the LSO crossing. In addition 
to the EOB resummed analytical expressions (dashed 
curves, red online), we also show the PN-expanded (1PN- 
accurate) absorbed fluxes as computed by Taylor and 
Poisson [2], (see also Eq. (13) of [3]). They are given by 

32 

-^ 1PN (^) = y^ 15 / 2 (l + 3.T), (25) 

-.F£ 1PN (aO = § ^ 17/2 - (26) 

When plotting these expressions we use two different PN 
representations of x: either x = (dashed-line, black 
online) consistently with the EOB waveform, or x^ = 
f2 2 / 3 (dash-dotted line, blue online). The two expressions 
differ only well below the LSO due to the violation of the 
Kepler constraint during the plunge. 

Following observations can be made in Fig. 5. First, 
the PN expanded expressions clearly underestimate the 
absorbed flux in the strong-field regime. This is expected 
due to the structure of the pf m in the circular case. It has 
been shown in Ref. [3] (Fig. 3) that at x = 1/7 0.14 
the lPN-accurate is more than a factor of two smaller 
than the corresponding /0 2 2 num computed from numerical 
data. 

Second, the EOB resummed expression (with the fitted 
coefficients cf m ) shows a very good consistency with the 
exact angular momentum flux computed from the waves. 
The fractional difference is « 1% at the beginning of 
the inspiral growing up to « 3% at the LSO crossing. 
Notably, this occurs also for the m = 1 flux, where the 



knowledge of the function comes completely from the 
fit to the circular data [3]. The fractional difference we 
find here is approximately one order of magnitude larger 
than for the asymptotic flux (for the same mass ratio v = 
10~ 3 ), sec Fig. 14 of [30]. This difference is not surprising 
because we have little analytical information to compute 
the EOB horizon flux. The computation relies mostly on 
the coefficients cf m obtained from the fit to the numerical 
data. 

Third, the fluxes stay close also below the LSO cross- 
ing, even though we do not expect the RWZ fluxes to be 
accurate close to the light-ring crossing. The fact that 
the fluxes remain so close during the late inspiral up to 
the plunge is by itself a confirmation that the fitted c,'s 
yield a rather accurate approximation to the coefficients 
one would get from the analytic PN calculation. 

In conclusion we have shown that the analytical , 
built using several pieces of information coming from a 
circularized binary (cither analytical or numerical) shows 
an excellent agreement with the exact horizon flux com- 
puted from the RWZ waves. This makes us confident 
that we can safely use as a new term in the radiation 
reaction to take horizon absorption into account. The 
influence of this term on the waveform phasing will be 
discussed in detail below. 



C. Effect on BBH phasing 

In this section we discuss and quantify the effect of 
the inclusion of absorbed fluxes, , in the dynamics 
on the observable GW (i.e. at infinity) from coalescing 
nonspinning binaries of different mass ratios. We work 
here only with EOB-gcnerated waveforms. 
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FIG. 6: (color online) Test-mass limit (y = 10 -3 ): including 
(v ; 0) in the dynamics and its effect on the £ = m = 2 
phasing. The top panel compares the I = m — 2 EOB wave- 
forms with (solid line) and without (dashed line) (v; 0). 
The accumulated phase difference (bottom panel) is of or- 
der 0.1 rad at LSO crossing (dash-dotted vertical line), and 
reaches a remarkable 1.5 rad at merger (dashed vertical line). 



We focus first on the test-mass limit, v — 10~ 3 , sub- 
ject to leading-order (in v) radiation reaction, Eq. (24) 
(we neglect then all the higher-order ^-dependent cor- 
rections). The effect of f H (v; 0) on the I = m = 2 
phasing is illustrated in Fig. 6. The initial separation 
is, as before, rg = 7, which yields about 41 orbits up to 
merger (see Table II). The top panel displays the EOB 
waveform without including horizon absorption (dashed 
line) together with the one where BH absorption is taken 
into account. The leftmost vertical line marks the LSO 
crossing, while the rightmost vertical line the light-ring 
crossing. The visible difference between the two wave- 
forms is made quantitative in the bottom panel of the 
figure, where the phase difference is shown. Here it is 



A022 = 4> H+j: — <j> . One sees that the phase difference 
is 0.1 rad at the LSO and grows up to 1.6 rad at merger. 

We turn now to compare a set of GWs from bina- 
ries with q = 1,4,10,50,100 and 1000, computed using 
the complete EOB dynamics. We run the simulations 
with and without horizon absorption and we compute 
the phase differences. The initial separation for q = 1000 
is ro = 7, while for the other mass ratios it is tq = 15, 
corresponding to initial GW frequency M<jj\ 2 = 0.0344. 
The result of this comparison is displayed in Fig. 7 and 
it is completed quantitatively by Table II. In the four 
panels of Fig. 7, the vertical lines mark, respectively from 
the left, the adiabatic LSO crossing and the EOB-defined 
light-ring crossing, i.e. the conventional location of the 
merger. First of all, we notice that even in the equal-mass 
case, where absorption effects are smallest and where the 
system has a limited number of cycles, one gets a dc- 
phasing of the order of 5 x 10~ 3 rad at the EOB merger. 
Remarkably, this value is comparable to (or just a lit- 
tle bit smaller than) the uncertainty on the phase of the 
most accurate numerical simulations of (equal-mass, non- 
spinning) coalescing black-hole binaries currently avail- 
able [36, 53, 54]. 

For higher mass ratios the cumulative effect of a larger 
horizon absorption (acting over more GW cycles) pro- 
duces larger and nonnegligible dephasings. As listed in 
Table II, mass ratios of q ~ 10 to 100 accumulate (re- 
spectively) a dephasing of A^f ~ 0.06 to 0.6 rad 
at LSO which increases of a factor 3 near the light 
ring, Acj^f ~ 0.22 to 2.2 rad. The last two columns 
in Table II list the dephasings obtained using the non- 
resummcd (lPN-accurate) radiation reaction. Interest- 
ingly, using such an expression of the absorbed flux yields 
dephasings that are up to 30% smaller (q = 100) at 
merger than the EOB prediction, underestimating the 
actual effect of absorption. 

Since horizon absorption effects on phasing are rela- 
tively large, especially for q > 50, they may be relevant 
in template modeling for large-mass-ratio binaries. In 
particular, we focus on IMR binaries made by stellar- 
mass compact object (SMCO) and an intermediate mass 
black-hole (IMBH), (M A ,M B ) ~ (1, 50 500)M S , that 
are candidate sources for Advanced LIGO [24], and for 
the Einstein Telescope (ET) [21]. We perform an indica- 
tive calculation of the faithfulness A [55] of an EOB tem- 
plate without absorption effects in describing a waveform 
with absorption effects. Given two (real) waveforms, say 
h\ (with horizon absorption) and h 2 (without horizon 
absorption) the faithfulness functional [55] (also denoted 
with the symbol T [56]) is defined as 



^4[/ii, hi] = max 



jhi,h 2 ) 

mump 



(27) 



where the maximization is done over a relative time r 
and phase shift a between the waveforms, and 



/>oo 

(h!,h 2 ) = m / df 

JO 



h(m(f) 

Sn(f) ' 



(28) 
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FIG. 7: (color online) Accumulated phase difference due to horizon absorption for different mass ratios q as obtained from 
EOB evolutions. The vertical lines mark the crossing of the EOB-defined LSO (leftmost line) and of the EOB-defined light-ring 
(rightmost line). For all binaries, the initial separation is ro = 15, corresponding to Mui%\ = 0.0344. 



defines the Wiener scalar product between the two sig- 
nals. Here, S n {f) is the one-sided power spectral density 
of the detector noise, h(f) the (complex) Fourier trans- 
form of the signal, and \\h\\ — (fr, /i) 1 / 2 the norm as- 
sociated to the Wiener scalar product. The mass ratios 
considered were q = 10, 50, 71.4286 and 100, correspond- 
ing to total masses M = (10 + 100)M o , (10 + 5OO)M , 
(1.4 + 1OO)M , and M = (14 + 14O)M . We followed 
the technical steps of Ref. [56] to compute accurately the 
Fourier transform of a EOB waveform. We computed 
A taking for S n both the ZER0_DET_HIGH_P anticipated 
sensitivity curve of Advanced LIGO [57] or that of the 
planned Einstein Telescope (ET) [58-60] . The numerical 
values of A are listed in Table III. Neglecting horizon 
absorption (for nonspinning binaries) leads to a loss of 
events (oc A 3 ) of, at most, 0.27% (for LIGO) and of, at 
most, 0.9% for ET. These numbers can be considered 
negligible for practical purposes. 



V. CONCLUSIONS 

We investigated the importance of horizon absorption 
effects in modelling GWs from nonspinning coalescing 
black hole binaries. Considering a recently proposed 
EOB resummed expression of the absorbed flux [3], we 
verified the EOB expression against pcrturbative wave- 
forms from large mass ratio (q = 1000) binaries (Sec. IV), 
and explored the effects of absorbed fluxes on the phas- 
ing considering EOB evolutions for binaries of different 
mass ratios q = 1 to 1000 (Sec. IV C). 

We tested the accuracy of the analytically resummed 
horizon flux [3], and in particular of the residual am- 
plitude corrections pf m , in the large- mass-ratio, pertur- 
bative limit. We compared it to the actual horizon 
flux of angular momentum computed solving the Regge- 
Wheeler-Zerilli equations in the time-domain. 

To improve the accuracy of the pcrturbative computa- 
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TABLE II: Accumulated phase differences due to to horizon absorption for different mass ratios. The data of the first six binaries 
are obtained from a complete EOB simulations. On the contrary, the dynamics of the last binary, shown for comparison, is that 
of a point-particle driven by leading-order radiation reaction only. For the first five binaries, the initial separation is ro = 15, 
which corresponds to frequency Mlo 2 , 2 w 0.0344, while the last two binaries start at ro = 7, i.e. Mlo 22 = 0.108. From left to 
right, the columns report: the mass ratio q; the symmetric mass ratio v — q/(l + q) 2 (y = 1/q for the last binary); the initial 
separation; the number of orbits up to merger (EOB-defined light-ring crossing), A or b; the dephasing A022 = </>22 +J * — ^22 
accumulated at the (adiabatic) EOB-defined LSO crossing; the corresponding value expressed in GW cycles; the dephasing 
accumulated at the EOB-defined light-ring crossing; the corresponding value expressed in GW cycles. The rightmost two 
columns show the phase difference accumulated using Taylor- Poisson, nonresummed, 1PN accurate radiation reaction. Note 
that the effect of horizon absorption on the phasing is still nonegligible (for q > 10) even using this leading order approximation 
to P$. 



q 


V 


ro 


A^orb 


A^f° [rad] 




A$ a a [rad] 


AJV LR 


A 1PN ^f° [rad] 


A 1PN ^ 2 R [rad] 


1 


0.250000 


15 


15 


0.003289 


0.000523 


0.005475 


0.000871 


0.002849 


0.004547 


4 


0.160000 


15 


21 


0.028725 


0.004572 


0.104712 


0.016665 


0.012320 


0.020246 


10 


0.082645 


15 


38 


0.064372 


0.010245 


0.220496 


0.035093 


0.052834 


0.199428 


50 


0.019223 


15 


153 


0.312210 


0.049690 


1.115319 


0.177508 


0.230220 


0.765105 


100 


0.009803 


15 


296 


0.620662 


0.098781 


2.217042 


0.352853 


0.458168 


1.549226 


1000 


0.000998 


7 


41.2 


0.129978 


0.020687 


1.453992 


0.231410 




1002 


0.000996 


7 


40.9 


0.129023 


0.020535 


1.563971 


0.248914 





TABLE III: Faithfulness between signals with and without 
horizon flux for SMCO-IMBH (nonspinning) binaries in the 
Advanced LIGO and Einsten Telescope sensitivity band. The 
merger frequency /merger corresponds to the maximum of the 
EOB waveform modulus | A22 1 



Q 


M a + M b [Mq] 


/"merger [Hz] 


-4aLIGO 


Aet 


10 


10 + 100 


89.16 


0.9999 


0.9998 


50 


10 + 500 


17.92 


0.9991 


0.9995 


71.43 


1.4 + 100 


89.21 


0.9991 


0.9983 


100 


1.4 + 140 


63.63 


0.9992 


0.9970 



tion, we employed two transmitting layers [31] (horizon- 
penetrating and hypcrboloidal) attached to a compact 
domain in standard Schwarzschild coordinates. This 
technique, summarized in Sec III, allows us to include 
in the computational domain both null-infinity, , and 
the horizon, H, via compactification in the tortoise coor- 
dinate. The resulting improvements of our perturbativc 
time-domain code combined with high-order finite differ- 
encing lead to such accurate computations of the inspiral 
and plunge that the late-time tail of the signal can be 
calculated very efficiently as reported in Appendix A. 

We computed the absorbed GW fluxes from the tran- 
sition from inspiral to plunge down to the late inspiral up 
to merger for the first time. We found that the quadrupo- 
lar contributions dominate over the subdominant multi- 
poles accouting for about 98% of the absorbed radiation 
(see bottom panel of Fig. 4). The £ = 2 absorbed an- 
gular momentum flux from the perturbative simulations 



proved to be consistent at the 1% level with the ana- 
lytical expressions proposed in [3]. Notably, the agree- 
ment remains excellent also below the LSO crossing and 
during the plunge. The resummation procedure for the 
flux introduced in [3] and the numerical determination 
of the higher-order PN terms entering the pf m amplitude 
corrections were crucial to obtain this result. The 1PN 
accurate, Taylor-expanded expression of the horizon flux 
as computed by Taylor and Poisson [2], underestimates 
horizon absorption by as much as a factor 2 during the 
late-inspiral and plunge phases. 

The absorbed flux of [3] has been used to build an 
additional term to the radiation reaction force of the 
EOB model, J 7 ^ , thereby incorporating in the model, 
in a resummed way, horizon absorption. By means of 
EOB simulations we explored its effect on the phasing of 
the GW emitted by binaries of different mass ratios q. 
Even in the current nonspinning case, it yields nonneg- 
ligiblc phase differences for q > 1. In particular, in the 
mass-ratio range q = 10 to 100 (see Table II), the accu- 
mulated phase differences are of the order 0.2 to 2 rad 
up to merger for circularized binaries initially at relative 
separation of ro = 15. By contrast, the PN-expandcd 
radiation reaction underestimates the dephasing by 9% 
to 48% (depending on q). 

Finally, we have performed a preliminary investigation 
of the impact of horizon absorption on the accurate mod- 
eling of templates for IMR nonspinning binaries made by 
a SMCO and a IMBH (M A , M B ) ~ (1, 50 + 500)M o . We 
found that neglecting would yield a loss of events by 
0.27% for Advanced LIGO and by 0.9% for ET. These 
losses are essentially negligible by current accuracy stan- 
dards. 

Horizon absorption effects are expected to be more rel- 
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evant for spinning binaries. It will be necessary to include 
them in J 7 ^ , after a suitable resummation procedure, to 
study their impact on the phasing. There is an ongoing 
study of this effect in the spinning case [17]. 
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Appendix A: Late-time tail decay for radial infall 
and insplunge trajectories 

In this Appendix we present, for the first time, the 
accurate computation of the late-time power-law tail of 
the waveform at generated by a particle plunging, 
both radially and following an inspiralling trajectory, into 
a Schwarzschild black hole. This result completes the 
knowledge of the ^-waveform for these events, already 
computed elsewhere [29, 30]. 

We recall that the gravitational waveform is computed 
by solving the RWZ equations in the time domain for 
each multipolc. The (5-function representing the particle 
is approximated by a narrow Gaussian of finite width 
a <C M, Eq. (Al). The representation of a particle 
as a Gaussian is a standard method when gravitational 
perturbations arc computed using finitc-diffcrcncc, time- 
domain methods. This representation, however, was con- 
sidered problematic, because time-domain codes gave rel- 
atively inaccurate results for gravitational fluxes [GI- 
GS]. Therefore, different prescriptions have been exper- 
imented with to improve on the representation of the 
point particle through a Gaussian [64-66]. Neverthe- 
less, the accuracy of time domain codes remained low, 
especially when compared with frequency domain ones. 
One open problem was the calculation of tail decay rates 
for a particle radially infalling into a Schwarzschild black 
hole [61]. 

Recently, a multi-domain hybrid method of finite dif- 
ference and spectral discretizations has been developed to 
solve this problem [67]. With this method, and using a 
large computational domain, the polynomially decaying 
part of the signal could be computed. However, the width 
of the Gaussian used in [67] to represent the particle is 
inadequate for the particle limit. In fact, the "particle" 
in this study is larger than the Schwarzschild black hole 
that provides the background. 



In this Appendix, we show that the accuracy provided 
by transmitting layers, combined with high-order finite 
differencing, allows us to calculate the tail decay rates 
accurately for realistic representations of a point particle 
in Schwarzschild spacetime. We present the decay rates 
not only for a radially infalling particle, but also for an 
insplunging one. 

As in previous work [26, 28-30], we approximate the 
delta distribution that represents the particle at time- 
dependent location, R*(t), by a Gaussian 

S(r ' - v^-P 1 ^)' (A1) 

Our prescription for the standard deviation, a, depends 
on resolution. We set a = 4 Ar* , so that the Gaussian is 
resolved well on our finite difference grid. 

Transmitting layers play an essential role in resolving 
narrow Gaussians because they allow us to compute the 
infinite domain solution in a small grid. This implies 
that the numerical resolution is not wasted in simulating 
empty space; instead, it can be focusscd to where the 
particle is located. As a consequence, we can afford to 
choose Ar„, and therefore the width of the Gaussian a, 
very small. 

Another advantage of using the layer method is that 
the implementation of high-order finite differencing be- 
comes simpler because there are no boundary conditions 
to be applied at either end of the domain. Note that 
even when good boundary conditions are available, their 
discretization and numerical implementation may not be 
straightforward. When no boundary conditions need to 
be applied, however, using a high order finite difference 
method becomes just a matter of widening the stencils. 

Using the transmitting layers, we have improved the 
accuracy of our previous work [30]. We use a smaller 
domain of [—20, 20] with interfaces at R± = ±12. Com- 
pared to our previous domain of [—50, 70], this gives us 
a factor of 3 in efficiency 7 . In addition, we use 8th order 
finite differencing as opposed to 4th order in [30]. As a 
result, we can compute the tail decay rates accurately, as 
reported below. 

1. Radial infall 

The calculation of gravitational perturbations caused 
by a particle falling radially into a non-rotating black hole 
is a classical problem in relativity [68, 69]. It serves as a 
good test bed for numerical computations, and there are 
still relatively recent studies on the problem [61, 67, 70]. 

We solve the radial infall of a particle to demonstrate 
the accuracy of our infrastructure. For a detailed de- 



7 By construction, reducing domain size does not decrease the time 
step for a given resolution. We did not attempt to find the opti- 
mal thickness for the layers. 
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FIG. 8: The evolution of the field for a radially infalling parti- 
cle starting at ro — 7. The plot spans 13 orders of magnitude. 
Observers are located (from top to bottom) at ,_0 , 30, and 15. 



scription of the setup, the reader is referred to the liter- 
ature [39, 50, 51, 71]. 

In Fig. 8, we show the absolute value of the Zerilli func- 
tion "020 caused by an infalling particle initially at rest at 
ro = as measured by three observers. The particle is 
represented by the Gaussian (Al) with a full width at half 
maximum (FWHM) 8 of 0.04M. We use 10,000 grid cells 
and a time stepping factor of 0.75 for the computation. 
Note that, differently from Refs. [50, 51] we put ■020 = 
initially and we do not solve consistently the Hamiltonian 
constraint. Since we are interested here in the late-time 
behavior of the waveform, this simplifying choice has no 
influence on our results. We see the QNM ringing after 
the plunge of the particle into the black hole, followed 
by late-time decay. The three curves in the figure corre- 
spond to the measurements of three observers (from top 
to bottom): the observer at infinity, the finite distance 
observer at 30M, and at 15M. The perturbations are 
computed for about 1000A/ which leads to a drop in the 
absolute value of the perturbation by 13 orders of mag- 
nitude. The polynomially decaying signal is reproduced 
accurately. 

The gain in accuracy is partly a result of the 8th order 
finite differencing, but mostly due to the high resolu- 
tion we can afford using transmitting layers, which allow 
us not only to compute the perturbations as measured 
by the observer at infinity, but also to follow the signal 
much longer than is possible with standard methods. For 
example, in Rcf. [67] the authors compute the perturba- 
tions until about 600A/ for a Gaussian source that has 
a FWHM of 5 — 10M which is larger than the size of 
the central black hole, and therefore cannot represent a 
realistic particle 9 . 



8 The FWHM of a normal distribution is given by its standard 
deviation a as 2<r\/2 In 2. 

9 The representation of the Gaussian in [67] leads to a FWHM of 
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FIG. 9: The local decay rates for the above evolution. The 
observers are located approximately at (from top to bottom) 
in units of M : { J?, 250, 80, 50, 35, 25, 20}. The dashed lines 
indicate the theoretically expected asymptotic decay rates: 
—4 at £ , and —7 at finite distances. 
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FIG. 10: The evolution of the real (solid line) and imaginary 
(dashed line) of the field for insplunge from ro = 7. The 
evolution spans 14 orders of magnitude. Observers are located 
(in units of M, from top to bottom) at J?, 35, and 18. 



We also plot the local decay rates as measured by dif- 
ferent far away observers in Fig. 9. The local decay rate 
plot gives a clear image of the accuracy of our compu- 
tation. We see that the expected decay rates are repro- 
duced accurately. The observer at infinity measures a 
rate of —4, whereas the rate for finite distance observers 
approaches —7. The intermediate behavior for the decay 
rates for these observers is in accordance with computa- 
tions of vacuum perturbations [41]. 

The local rates for the observers at 25M and 20M in 
Fig. 9 have been cut from the plot at late times because of 
large oscillations. The loss of accuracy for these observers 



2\J a In 2. The authors present studies with a ranging between 
10 and 50. 
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is not only because of accumulated truncation error, but 
mostly because the fast decaying signal reaches machine 
precision. If necessary, the decay rate calculation can 
be further improved by using quadruple precision, and 
possibly higher resolution. 

2. Insplunge 

The main interest in this paper is the study of particles 
plunging into the central black hole following a phase of 
quasi-circular inspiral (insplunge). We compute the tail 
decay rates also for this case. As above, the initial sep- 



aration is ro = 7. In Fig. 10 we show the absolute value 
of the real part (solid line) and imaginary part (dashed 
line) of the perturbation, again as measured by three ob- 
servers (from top to bottom): the observer at infinity and 
the finite distance observers at 35M and 18M. The com- 
putational parameters are the same as in the radial infall 
study. We see that the field is followed for 14 orders of 
magnitude, and the evolution is presented until 1500A/ 
this time. The three stages of the evolution (inspiral, 
ringing, and polynomial decay) are clearly visible. The 
local decay rates show qualitatively the same behavior as 
in Fig. 9 and are therefore not plotted. 
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